{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 16\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code from previous notebooks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def update_func(state, t, system):\n",
    "    \"\"\"Update the thermal transfer model.\n",
    "    \n",
    "    state: State (temp)\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: State (temp)\n",
    "    \"\"\"\n",
    "    r, T_env, dt = system.r, system.T_env, system.dt\n",
    "    \n",
    "    T = state.T\n",
    "    T += -r * (T - T_env) * dt\n",
    "    \n",
    "    return State(T=T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_simulation(system, update_func):\n",
    "    \"\"\"Runs a simulation of the system.\n",
    "    \n",
    "    Add a TimeFrame to the System: results\n",
    "    \n",
    "    system: System object\n",
    "    update_func: function that updates state\n",
    "    \"\"\"\n",
    "    init = system.init\n",
    "    t_0, t_end, dt = system.t_0, system.t_end, system.dt\n",
    "    \n",
    "    frame = TimeFrame(columns=init.index)\n",
    "    frame.row[t_0] = init\n",
    "    ts = linrange(t_0, t_end, dt)\n",
    "    \n",
    "    for t in ts:\n",
    "        frame.row[t+dt] = update_func(frame.row[t], t, system)\n",
    "    \n",
    "    return frame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(T_init, r, volume, t_end):\n",
    "    \"\"\"Makes a System object with the given parameters.\n",
    "\n",
    "    T_init: initial temperature in degC\n",
    "    r: heat transfer rate, in 1/min\n",
    "    volume: volume of liquid in mL\n",
    "    t_end: end time of simulation\n",
    "    \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    init = State(T=T_init)\n",
    "                   \n",
    "    return System(init=init,\n",
    "                  r=r, \n",
    "                  volume=volume,\n",
    "                  temp=T_init,\n",
    "                  t_0=0, \n",
    "                  t_end=t_end, \n",
    "                  dt=1,\n",
    "                  T_env=22)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `root_bisect`\n",
    "\n",
    "As a simple example, let's find the roots of this function; that is, the values of `x` that make the result 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def func(x):\n",
    "    return (x-1) * (x-2) * (x-3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`modsim.py` provides `root_bisect`, which searches for a root by bisection.  The first argument is the function whose roots we want.  The second argument is an interval that contains a root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = root_bisect(func, [0.5, 1.5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result is an object that contains the root that was found and other information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we provide a different interval, we find a different root."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = root_bisect(func, [1.5, 2.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the interval doesn't contain a root, the results explain the error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = root_bisect(func, [4, 5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to find the value of `r` that makes the final temperature 70, so we define an \"error function\" that takes `r` as a parameter and returns the difference between the final temperature and the goal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def error_func1(r):\n",
    "    \"\"\"Runs a simulation and returns the `error`.\n",
    "    \n",
    "    r: heat transfer rate, in 1/min\n",
    "    \n",
    "    returns: difference between final temp and 70 C\n",
    "    \"\"\"\n",
    "    system = make_system(T_init=90, r=r, volume=300, t_end=30)\n",
    "    results = run_simulation(system, update_func)\n",
    "    T_final = get_last_value(results.T)\n",
    "    return T_final - 70"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `r=0.01`, we end up a little too warm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_func1(r=0.01)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With `r=0.02`, we end up too cold."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "error_func1(r=0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The return value from `root_bisect` is an array with a single element, the estimated value of `r`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = root_bisect(error_func1, [0.01, 0.02])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_coffee = res.root"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we run the simulation with the estimated value of `r`, the final temperature is 70 C, as expected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)\n",
    "results = run_simulation(coffee, update_func)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:**  When you call `root_bisect`, it calls `error_func1` several times.  To see how this works, add a print statement to `error_func1` and run `root_bisect` again."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Repeat this process to estimate `r_milk`, given that it starts at 5 C and reaches 20 C after 15 minutes.  \n",
    "\n",
    "Before you use `root_bisect`, you might want to try a few values for `r_milk` and see how close you can get by trial and error.  Here's an initial guess to get you started:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_milk = 0.1\n",
    "milk = make_system(T_init=5, r=r_milk, volume=50, t_end=15)\n",
    "results = run_simulation(milk, update_func)\n",
    "T_final = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Mixing liquids"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function takes `System` objects that represent two liquids, computes the temperature of the mixture, and returns a new `System` object that represents the mixture."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mix(system1, system2):\n",
    "    \"\"\"Simulates the mixture of two liquids.\n",
    "    \n",
    "    system1: System representing coffee\n",
    "    system2: System representing milk\n",
    "    \n",
    "    returns: System representing the mixture\n",
    "    \"\"\"\n",
    "    assert system1.t_end == system2.t_end\n",
    "    \n",
    "    V1, V2 = system1.volume, system2.volume\n",
    "    T1, T2 = system1.temp, system2.temp\n",
    "    \n",
    "    V_mix = V1 + V2\n",
    "    T_mix = (V1 * T1 + V2 * T2) / V_mix\n",
    "    \n",
    "    return make_system(T_init=T_mix,\n",
    "                       r=system1.r,\n",
    "                       volume=V_mix,\n",
    "                       t_end=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`mix` requires the `System` objects to have `temp` as a system variable.  `make_system` initializes this variable;\n",
    "the following function makes sure it gets updated when we run a simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_and_set(system):\n",
    "    \"\"\"Run a simulation and set the final temperature.\n",
    "    \n",
    "    system: System\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    results = run_simulation(system, update_func)\n",
    "    system.temp = get_last_value(results.T)\n",
    "    return results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixing immediately\n",
    "\n",
    "Next here's what we get if we add the milk immediately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)\n",
    "coffee.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "milk = make_system(T_init=5, r=r_milk, volume=50, t_end=30)\n",
    "milk.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_first = mix(coffee, milk)\n",
    "mix_first.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_results = run_and_set(mix_first)\n",
    "mix_first.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixing at the end\n",
    "\n",
    "First we'll see what happens if we add the milk at the end.  We'll simulate the coffee and the milk separately."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "coffee_results = run_and_set(coffee)\n",
    "coffee.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "milk_results = run_and_set(milk)\n",
    "milk.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the results look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(coffee_results.T, label='coffee')\n",
    "plot(milk_results.T, '--', label='milk')\n",
    "\n",
    "decorate(xlabel='Time (minutes)',\n",
    "         ylabel='Temperature (C)',\n",
    "         loc='center left')\n",
    "\n",
    "savefig('figs/chap16-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what happens when we mix them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_last = mix(coffee, milk)\n",
    "mix_last.temp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "mix_last.temp - mix_first.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following function takes `t_add`, which is the time when the milk is added, and returns the final temperature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_and_mix(t_add, t_total):\n",
    "    \"\"\"Simulates two liquids and them mixes them at t_add.\n",
    "    \n",
    "    t_add: time in minutes\n",
    "    t_total: total time to simulate, min\n",
    "    \n",
    "    returns: final temperature\n",
    "    \"\"\"\n",
    "    coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=t_add)\n",
    "    coffee_results = run_and_set(coffee)\n",
    "    \n",
    "    milk = make_system(T_init=5, r=r_milk, volume=50, t_end=t_add)\n",
    "    milk_results = run_and_set(milk)\n",
    "    \n",
    "    mixture = mix(coffee, milk)\n",
    "    mixture.t_end = t_total - t_add\n",
    "    results = run_and_set(mixture)\n",
    "\n",
    "    return mixture.temp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try it out with a few values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_and_mix(t_add=0, t_total=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_and_mix(t_add=15, t_total=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "run_and_mix(t_add=30, t_total=30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And then sweep a range of values for `t_add`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "sweep = SweepSeries()\n",
    "for t_add in linspace(0, 30, 11):\n",
    "    sweep[t_add] = run_and_mix(t_add, 30)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the result looks like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot(sweep, label='final temp', color='C2')\n",
    "decorate(xlabel='Time added (min)',\n",
    "         ylabel='Final temperature (C)')\n",
    "\n",
    "savefig('figs/chap16-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the analytic result to compute temperature as a function of time.  The following function is similar to `run_simulation`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_analysis(system):\n",
    "    \"\"\"Computes temperature using the analytic solution.\n",
    "        \n",
    "    system: System object\n",
    "    \n",
    "    returns: TimeFrame\n",
    "    \"\"\"\n",
    "    T_env, r = system.T_env, system.r\n",
    "    \n",
    "    T_init = system.init.T    \n",
    "    ts = linspace(0, system.t_end)\n",
    "    \n",
    "    T_array = T_env + (T_init - T_env) * exp(-r * ts)\n",
    "    \n",
    "    # to be consistent with run_simulation,\n",
    "    # we put the array into a TimeFrame\n",
    "    return TimeFrame(T_array, index=ts, columns=['T'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we run it.  From the analysis (see `chap16sympy.ipynb`), we have the computed value of `r_coffee2`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_coffee2 = 0.011610223142273859\n",
    "coffee2 = make_system(T_init=90, r=r_coffee2, volume=300, t_end=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "results = run_analysis(coffee2)\n",
    "T_final_analysis = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can compare to the results from simulation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "coffee = make_system(T_init=90, r=r_coffee, volume=300, t_end=30)\n",
    "results = run_simulation(coffee, update_func)\n",
    "T_final_simulation = get_last_value(results.T)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They are identical except for a small roundoff error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_final_analysis - T_final_simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "**Exercise:**  Suppose the coffee shop won't let me take milk in a separate container, but I keep a bottle of milk in the refrigerator at my office.  In that case is it better to add the milk at the coffee shop, or wait until I get to the office?\n",
    "\n",
    "Hint: Think about the simplest way to represent the behavior of a refrigerator in this model.  The change you make to test this variation of the problem should be very small!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
